Importance of cooling in triggering the collapse of hypermassive neutron stars 
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The inspiral and merger of a binary neutron star (NSNS) can lead to the formation of a hy- 
permassive neutron star (HMNS). As the HMNS loses thermal pressure due to neutrino cooling 
and/or centrifugal support due to gravitational wave (GW) emission, and/or magnetic breaking of 
differential rotation it will collapse to a black hole. To assess the importance of shock-induced ther- 
mal pressure and cooling, we adopt an idealized equation of state and perform NSNS simulations 
in full GR through late inspiral, merger, and HMNS formation, accounting for cooling. We show 
that thermal pressure contributes significantly to the support of the HMNS against collapse and 
that thermal cooling accelerates its "delayed" collapse. Our simulations demonstrate explicitly that 
cooling can induce the catastrophic collapse of a hot hypermassive neutron star formed following 
the merger of binary neutron stars. Thus, cooling physics is important to include in NSNS merger 
calculations to accurately determine the lifetime of the HMNS remnant and to extract information 
about the NS equation of state, cooling mechanisms, bar instabilities and B-fields from the GWs 
emitted during the transient phase prior to BH formation. 

PACS numbers: 04.25.D-,04.25.dk,04.30.-w 



I. INTRODUCTION 

The inspiral and merger of compact binaries has at- 
tracted considerable attention in recent years for two 
main reasons. First, such systems emit a large flux 
of gravitational waves (GWs), making them among the 
most promising sources for GWs detectable by ground- 
based laser interferometers such as LIGO 0, Hjj, VIRGO 
0,13, GEO @, and KAGRA @, as well as by proposed 
space-based interferometers such as eLISA/NGO [7j and 
DECIGO 0. Second, black hole - neutron star (BHNS) 
and neutron star - neutron star (NSNS) mergers are can- 
didates for the central engines that power the observed 
short-hard gamma ray bursts (sGRBs). 

Extracting physical information about these binaries 
from their GWs and their accompanying electromagnetic 
signals may reveal critical details about the equation of 
state of neutron star matter and may unveil the na- 
ture of the sGRB phenomenon. However, interpreting 
the data requires careful modeling of these systems in 
full general relativity (see Q for a comprehensive review 
and references). Most effort in general relativity to date 
has focused on modeling black hole-black hole (BHBH) 
binaries (see also [I3|), and neutron star-neutron star 
(NSNS) binaries (see also [ll|), with some recent work 
on black hole-neutron star binari es ( see also [HJ), and 
white dwarf-neutron star binaries [13h15| . 

NSNSs are known to exist, which makes NSNS systems 
particularly attractive to study. Theoretical calculations 
show that NSNS mergers can lead to the formation of 
a hypermassive neutron star. A HMNS 16] is a differ- 
entially rotating NS whose mass exceeds the maximum 
mass of a uniformly rotating star [17|, [l8| . The latter is 
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about 20% larger than the maximum mass of a nonro- 
tating (spherical) equilibrium star (the TOV limit) [Tfj| . 
Typically a HMNS forms following the merger of a NSNS, 
when the system's total mass is smaller than some thresh- 
old mass Mth- According to [19j] this threshold mass is 
M th « 1.3 - 1.35M sph , where M sph is the TOV limit for 
the same EOS. 

A HMNS is a transient, quasicquilibrium configura- 
tion. It will eventually undergo "delayed collapse" on 
a secular (dissipative) time scale, which may power a 
sGRB. There are two distinct routes by which this col- 
lapse might be triggered: 

1. If the HMNS is primarily centrifugally supported, 
redistribution of an gula r momentum by viscosity or 
magnetic fields [2fJ, [21| , and /or loss of angular mo- 
mentum by GW emission [22j destroys the support 
provided, leading to catastrophic collapse. 

2. If the HMNS is primarily supported by thermal 
pressure generated by shocks during merger, de- 
layed collapse may be triggered by the loss via neu- 
trino cooling of thermal energy [231 ] . 

While catastrophic collapse of a cold HMNS via viscos- 
ity or magnetic fields has been demonstrated using fully 
general relativistic calculations [20j, |2l| , there are no fully 
general relativistic calculations to date that demonstrate 
explicitly that cooling can induce collapse of a hot HMNS 
produced following the merger of binary neutron stars. 

HMNSs formed in NSNS mergers will always be hot 
due to shock heating. A priori it is not clear which mech- 
anism is most important for holding up a HMNS against 
collapse: centrifugal forces or thermal pressure. The an- 
swer to this question is still open and may depend on the 
nature of the companions (e.g. masses, EOS etc.). 

Recent simulations of binary NS mergers that form hy- 
permassive NSs seem to point in different directions. For 
example, in [24|, HH an equal-mass NSNS is evolved as- 
suming af = 2 equation of state (EOS). It is shown that 
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angular momentum carried away by gravitational waves 
alone can induce the collapse. Reference [26[ also evolves 
an equal-mass NSNS, but with a more realistic, finite 
temperature, nuclear EOS. They find that the deviation 
of their HMNSs from axisymmetry is so small that GW 
emission is significantly reduced. The authors argue that 
shock heating is sufficiently important that their HMNSs 
are supported by the excess thermal pressure. 

Determining which mechanism controls the lifetime of 
the remnant is important because it determines the time 
interval between the NSNS merger and the delayed col- 
lapse - a time interval that can in principle be measured 
by Advanced LIGO/VIRGO. It is the time interval be- 
tween the end of the gravitational wave signal due to the 
inspiral and the beginning of the burst signal due to the 
delayed collapse. If differential rotation support is most 
important, then the time interval is governed by, e.g. the 
Alfven time scale, assuming magnetic braking of differ- 
ential rotation is most important, or the GW time scale, 
in the case of a rapidly spinning remnant that develops a 
bar. By contrast, if thermal pressure is dominant, then 
the time scale is governed by thermal cooling. Therefore, 
knowing the mechanism driving collapse may place con- 
straints on seed magnetic field magnitudes, or the exis- 
tence of bar modes, or the relevant cooling mechanisms. 
It could even place constraints on the temperature of 
matter, as well as the nuclear EOS. 

To disentangle the effects of thermal support from 
those of rotational support, previous studies compared 
results from NSNS simulations that suppress shocks (by 
enforcing a strictly cold EOS) to those that allow shocks. 
If the HMNS remnant lives longer with shocks than with- 
out, then it is tempting to infer that thermal pressure due 
to shock heating is chiefly responsible for supporting the 
remnant. However it is not possible to draw such a firm 
conclusion because shocks, which act on a hydrodynam- 
ical time scale, not only heat the gas, thereby increasing 
the total pressure support, but also affect the matter and 
angular momentum profiles. Different profiles can them- 
selves increase the lifetime of a HMNS. 

The goal of this paper is to study the relative impor- 
tance of thermal pressure in supporting HMNSs from col- 
lapse and demonstrate that cooling can induce the catas- 
trophic collapse of a HMNS formed following the merger 
of binary neutron stars. We accomplish this by perform- 
ing a limited set of NSNS simulations in full GR through 
late inspiral, merger, (hot) HMNS formation, and col- 
lapse. We account for cooling in the HMNS remnant via 
a covariant cooling scheme we developed in [f4| • We then 
compare this HMNS evolution to a control simulation, in 
which the cooling mechanism is disabled. 

Our simulations model the initial NSNS binary as 
equal-mass, irrotational, quasiequilibrium n = 1 poly- 
tropes in a quasicircular orbit, corresponding to case 
1. 46-45-* of 0- 

Following the NSNS merger, a quasiequilibrium HMNS 
forms. We then continue the evolution of the remnant 
with and without cooling, which we model via an effective 



local emissivity. For the runs with cooling we choose two 
cooling time scales. We find that, independent of the 
cooling time scale chosen, the HMNS collapses and forms 
a BH within a few cooling time scales. 

Our simulations suggest that shock-induced thermal 
pressure is a significant source of support against gravi- 
tational collapse, even in the case of polytropic NSs and 
demonstrate explicitly that cooling can induce the catas- 
trophic collapse of a HMNS. Estimating the temperature 
of the remnant, we find that a realistic neutrino cool- 
ing time scale is of order a few 1 00ms. Given that our 
estimated cooling time scale is comparable to the angu- 
lar momentum redistribution/loss time scales due to ei- 
ther magnetic braking or GWs, our results suggest that 
accounting for cooling is a critical ingredient in predict- 
ing the lifetime of a HMNS. Accordingly, cooling physics 
must be incorporated in models of binary NS simulations. 

The paper is structured as follows. In Sec. |H] we re- 
view the time scales relevant to HMNSs formed in binary 
NSNS mergers. Sections Mil and IIVI summarize the ini- 
tial data, basic evolution equations, numerical methods, 
and cooling formalism. The basic results are presented 
in Sec. |V] and summarized in Sec. I VII Throughout this 
work, geometrized units are adopted, where G = c = I , 
unless otherwise specified. 



II. TIME SCALES 

The relevant time scales in the evolution of a typical 
HMNS formed in NSNS mergers are its rotation period 
T, the gravitational wave time scale £gw, the cooling 
time scale t coo \, and Alfven time scale £a- We provide 
rough estimates of these time scales in this section. 



A. Rotation period 

We express the HMNS angular frequency f2 as some 
fraction e of the break-up angular frequency f2 ms 



(1) 



where M is the HMNS mass and R its radius. The rota- 
tion period of the HMNS can then be written as 
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For the numerical estimate we have used the values for 
the mass and radius of a typical HMNS remnant. 
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Gravitational wave time scale 



© in Eq. © we find 



GW emission sets the time scale of angular momentum 
loss from the system. The gravitational wave time scale 
for a triaxial, incompressible, spinning ellipsoid with el- 
lipticity e can be estimated as [271 ] 
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where J ps MR 2 Q, is the HMNS angular momentum and 
the ellipticity is defined as 
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where a is the semi- major axis of the HMNS, b the semi- 
minor axis, and R is (a+b)/2. To estimate the time scale, 
we assumed a value for the ellipticity that corresponds 
to a plausible bar. Note also that our estimated £qw 
is comparable to the GW time scale inferred by direct 
numerical simulations in 12511. 
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where n = p/m n , with p = 3M/4ttR 3 the mean HMNS 
density, and m n the mass of a neutron. For the numer- 
ical estimates above we used typical rms values for the 
neutrino energy of order lOMeV, as found in the simu- 
lations of 28]. Note that for typical neutrino energies 
of 20MeV found in [2(| the neutrino cooling timescale is 
~ 2s. Both of these works used approximate neutrino 
transfer schemes. We see that obtaining a neutrino cool- 
ing time scale depends on identifying the energy(ies) of 
typical neutrino (s), which in turn requires accurate mod- 
eling of not only bulk motion but also the microphysics. 



D. Alfven time scale 

Magnetic fields set the time scale for the braking of 
differential rotation in typical HMNSs. This occurs on 
the Alfven time scale [21| . given by 



C. Cooling time scale 



HMNSs are cooled predominantly by emission of neu- 
trinos. At densities > 10 11 g/cm 3 neutrinos become 
trapped [l?) • Therefore, the cooling time scale is set by 
the time it takes for the neutrinos to diffuse out of the 
hot HMNS remnant. The main sources of opacity are 
free nucleon scattering and neutrino absorption by nu- 
cleons (since protons and neutrons comprise the bulk of 
the HMNS). The diffusion time scale can be estimated as 



tc 



I*- 



(5) 



where A n is the mean free path of the neutrinos given by 

(6) 



A„ 1 = na n 



where n is the neutron number density [29j. a n is the to- 
tal interaction cross section er„ = <7 sca t + <r a bs , where the 
elastic scattering and absorption cross sections are re- 
spectively given by [13, HH 
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where <7q = 1-76 x 10 44 cm 2 , m e is the electron mass, 
and E v the neutrino energy. Substituting Eqs. ^ and 
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where va is the Alfven velocity, and where a strong but 
dynamically unimportant interior magnetic field has been 
assumed for the numerical estimate. While little is known 
about the strength of NS interior magnetic fields, the 
value appearing in ([9]) is consistent with magnetars mod- 
els [13] ■ In addition, NSNS simulations indicate that 
magnetic instabilities can amplify interior B-fields from 
- 10 12 G to - 10 15 G during merger [H|. 



E. Time scale summary 

These time scale estimates indicate that the neutrino 
cooling time scale can be comparable to the magnetic 
braking/ angular momentum loss time scales in typical 
HMNSs. If thermal pressure is the dominant source of 
support in an HMNS against catastrophic collapse to a 
BH, then the cooling time scale will determine the time 
interval between the GW signals at merger and collapse. 
Even if thermal pressure contributes only partially to the 
support of the HMNS, the remnant will collapse faster 
with cooling than without. These considerations neces- 
sitate the modeling of neutrino cooling in simulations 
of NSNS mergers that form HMNSs, not only to pre- 
dict the neutrino signature, but also to determine what 
mechanism drives the remnant to its final configuration. 
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Knowing the results from such simulations, it may be 
possible to extract useful information about the temper- 
ature of the matter, neutrino cooling mechanisms, the 
existence of bar modes, and the magnetic field strength 
and possibly place constraints on the nuclear EOS from 
the GW observations. We perform preliminary simula- 
tions to probe this issue below. 



III. BASIC EQUATIONS 

This section introduces our notation, summarizes our 
methods and numerical techniques as described in (32i — 
[35j . Greek indices denote all four spacetime dimensions 
(0, 1, 2, and 3), and Latin indices label spatial parts only 
(1, 2, and 3). 

We use the 3+1 formulation of general relativity and 
decompose the metric into the following form: 

ds 2 = -a 2 dt 2 + J ij (dx i + P i dt){dx j + frdt) . (10) 

The fundamental variables for metric evolution are the 
spatial three-metric 7^ and extrinsic curvature Kij. 
We adopt the Baumgarte-Shapiro-Shibata-Nakamura 
(BSSN) formalism [3a . |37| in which the evolution vari- 
ables are the conformal exponent <j> = ln(7)/12, the con- 
formal 3-metric 7^ = e~ 4 ^jij, three auxiliary functions 
V = — 7 lJ 1 j, the trace of the extrinsic curvature K, and 
the trace-free part of the conformal extrinsic curvature 
A i:j = e-^>(Kij - jijK/3). Here, 7 = det( 7ij ). The full 
spacetime metric g^ v is related to the three-metric 7 M „ 
by 7(ny = g^v + n IJt n v , where the future-directed, timclikc 
unit vector ti} 1 normal to the time slice can be written in 
terms of the lapse a and shift /3 l as tt/ 1 = oT l (1, — f3 l ). 
Evolution equations for these BSSN variables are given 
by Eqs. (9)-(13) in [32j. We adopt the standard punc- 
ture gauge conditions: an advective "1+log" slicing con- 
dition for the lapse and a 'T-freezing" condition for the 
shift [38j]. The evolution equations for a and fP are 
given by Eqs. (2)-(4) in [33[, with the 77 parameter set 
to 0.2/M, where M is the ADM mass of the NSNS bi- 
nary. We add a fifth-order Kreiss-Oliger dissipation term 
to all evolved BSSN, lapse and shift variables to reduce 
high-frequency numerical noise associated with AMR re- 
finement interfaces. 

The fundamental hydrodynamic (HD) variables are the 
rest-mass density po, specific internal energy e, pressure 
P, and four- velocity it M . We adopt a T-law equation of 
state (EOS) P = (T - l)p e with r = 2, which reduces 

to an n = 1 polytropic law [P = Kp^ +1 ^ n ^] for the initial 
(cold) neutron star matter. The fluid stress-energy tensor 
is given by 

T^lv = Pvhu^Uv + Pg^ u , (11) 

where h = 1 + e + P/po is the specific enthalpy. 

In the standard numerical implementation of the gen- 
eral relativistic hydrodynamic (GRHD) equations using 



a conservative scheme, it is useful to introduce the "con- 
servative" variables p*, Si, f. They are defined as 

P* = -y/lPon^u' 1 , (12) 

§i = -VrWy; , (13) 

f = ^T^n v - p* . (14) 

The evolution equations for /?*, Si and f can be derived 
from the conservation of rest mass V At (p*w M ) = and the 
conservation of energy-momentum V^T^ = 0, giving 
rise to Eqs. (27)-(30) in 0. 



IV. NUMERICAL METHODS 

A. Initial data 

For initial data we choose an irrotational NSNS sys- 
tem in a quasiequilibrium circular orbit that consists of 
equal-mass, n — 1 polytropic NSs. The initial data sat- 
isfy the conformal thin sandwich equations 9] , have been 
calculated using the LORENE spectral methods numerical 
libraries [39( and are publicly available. These data apply 
to a configuration with arbitrary k, compaction (in iso- 
lation) M/R = 0.12, where the compaction of the max- 
imum mass configuration is M/R = 0.216. Each star 
has a rest mass that is 72% of the maximum allowable 
TOV rest mass for this EOS. The initial cold configura- 
tion has a coordinate separation of 11.31M, where M is 
the ADM mass of system, with Mil — 0.024, where is 
the angular frequency of the system. The ADM angular 
momentum of the system is J/M 2 = 1.02. We note here 
that our initial data correspond to case 1. 46-45-* of [24| 
and that they can be considered as the polytropic coun- 
terpart of case H studied in (26|. If we set k — 393.9 km 2 , 
the ADM mass of our stars in isolation becomes 1.59M , 
which is very close to the ADM mass (1.6M Q ) in isola- 
tion of case H in (26j, where a finite temperature EOS 
was adopted that yields for zero-temperature matter a 
maximum TOV mass of 2.2M Q . 



B. Evolution of the metric and matter 

We evolve the BSSN equations with fourth-order ac- 
curate, centered finite-differencing stencils, except on 
shift advection terms, where we use fourth-order accurate 
upwind stencils. We apply Sommcrfeld outgoing wave 
boundary conditions to all BSSN fields. Our code is em- 
bedded in the CACTUS parallelization framework [40| . 
and our fourth-order Runge-Kutta timcstepping is man- 
aged by the MoL (Method of Lines) thorn, with a Courant- 
Friedrichs-Lewy (CFL) factor set to 0.45. We use the 
Carpet [4l| infrastructure to implement the moving-box 
adaptive mesh refinement. In all AMR simulations pre- 
sented here, we use second-order temporal prolongation, 
coupled with fifth-order spatial prolongation. 
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FIG. 1. Case A orbital-plane rest-mass density contours at selected times. Contours are plotted according to po = 
po,max(10 _0,37j '' _0 131 ), (j=0, 1, ... 8). The color sequence dark red, red, orange, yellow, green, light green, blue and 
light blue implies a sequence from higher to lower values. This roughly corresponds to darker grey-scaling for higher 
values. The maximum initial NS density is Atpo,max = 0.0917, or po.max = 4.58 x 10 g cm - 3 (1.45M /M O ) 2 . Here 
M = 1.32 x 10- 5 (M /1.45A'/ Q )s= 3.98(M /1.45M a )km is the ADM mass, and Mq denotes the rest mass of each star. 



The GRHD equations are evolved via a high-resolution 
shock-capturing (HRSC) technique 42] that employs 
PPM [43[ coupled to the Harten, Lax, and van Leer 
(HLL) approximate Riemann solver 01 . The adopted 
GRHD scheme is second-order accurate for smooth flows, 
and first-order accurate when discontinuities (e.g. shocks) 
arise. To stabilize our scheme in regions where there is 
no matter, we maintain a tenuous atmosphere on our 
grid, with a density floor p atm set equal to 10~ 10 times 
the initial maximum density on our grid. The initial 
atmospheric pressure P a t m is set equal to the cold poly- 
tropic value Patm = K Patm- Throughout the evolution, 
we impose limits on the atmospheric pressure to pre- 
vent spurious heating and negative values of the inter- 
nal energy e due to numerical errors. Specifically, we 
require P min < P < P max , where P max = IQnp^ and 
Pnin = K /°o/2- Whenever P exceeds P max or drops be- 
low P mm , we reset P to P max or P m i n , respectively. We 
impose these pressure limits only in regions where the 
rest-mass density remains very low (po < 100p a t m ), as 
in (al. 



C. Radiative cooling 

We now briefly describe our method for implementing 
cooling in our simulations. For a derivation and details 
regarding this covariant cooling method, see Q3 ■ 

The dynamics of radiation is governed by [45l447j | 

V a R al3 = -G p , (15) 
where R a/3 is the radiation stress-energy tensor given by 

R af3 = J dvdVLl v N a N^ (16) 

and G a is the radiation four-force density given by 

G a — f dudn( X uh-j v )N a . (17) 



In the equations above dVL is the solid angle, v and 
Iv = Iv{x a , N l , v) are the radiation frequency and spe- 
cific intensity of radiation at x a moving in direction 
N a = p a /hv, respectively. All quantities are measured 
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FIG. 2. Case A meridional (XZ) plane rest-mass density 
(upper panel) and K contours (lower panel). Density con- 
tours are plotted according to po = po,max(10 -0 ' 375 -' -0 ' 131 ), 
(i=0, 1, ... 8), where Kpo,max = 0.0917, or p ,max = 4.58 x 
10 14 g cm _3 (1.45M0/Mo) 2 . K contours are plotted according 
to K = Jf max 10- 031j , 0=0, 1, ... 8). Here A~ max = 1.77. 
The color coding is the same as used in Fig. [T] In the lower 
panel light blue indicates K w 1 and dark red K ~ 1.6. Here 
M = 1.32 x 10~ 5 (M /1.45M©)s= 3.98(M /1.45M o )km is the 
ADM mass, and Mq denotes the rest mass of each star. 



in the local Lorentz frame of a fiducial observer with four- 
velocity uf id , i.e., 

hv = -p a uf id , (18) 

where p a is the photon four-momentum and h denotes 
Planck's constant. The energy- momentum conservation 
equation then becomes 

V a {T aP + R a0 ) = 



(19) 



or after using Eq. ([15 



V a T al3 = G . 



(20) 



Our artificial cooling prescription amounts to finding 
a functional form for C 3 such that thermal energy and 
pressure are drained from the system. Choosing 



and setting 



G a = -u a A, 



A = — £th, 



(21) 



(22) 



where r c is some prescribed cooling time scale, it can be 
shown that in a frame comoving with the fluid the specific 
thermal energy of a fluid parcel evolves as follows [lij 
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(23) 



where r is the proper time of a comoving observer. 

The first term in brackets on the RHS of Eq. (|2"3")l arises 
from adiabatic compression or expansion. The second 
term corresponds to cooling and radiates away thermal 
energy exponentially. 

Projecting Eq. (|2TJ)) using the timelike unit vector n a 
normal to spacelike hypersurfaces and the projection op- 
erator h a fs = 5 a f} + n a ni3 , we find that the 3+1 GRHD 
equations become 

dtS, + djia^jT'i) = \^T ali g aSi ^ - a^A, (24) 
and 

dtf + ft(a 2 V7T 04 - P *v l ) = s - a 2 V7«°A, (25) 

where we have used Eq. (|21l) . Thus cooling enters as a 
source term in the GRHD equations. 



D. Recovery of primitive variables 

At each timestep, we need to recover the "primitive 
variables" po, P, and v l from the "conservative" variables 
p*, f, and Si. We perform the inversion by numerically 
solving two nonlinear equations via the Newton- Raphson 
method as described in [48j . using the code developed 
in HI. 

Sometimes the "conservative" variables may assume 
values which are out of physical range, resulting in un- 
physical primitive variables after inversion (e.g. negative 
pressure or even complex solutions). This usually hap- 
pens in the low-density "atmosphere" or deep inside the 
BH interior (when a BH is present) where high-accuracy 
evolution is difficult to maintain. Various techniques 
have been suggested to handle the inversion failure (see, 
e.g. Our approach is mainly to impose constraints 

on the conservative variables to reduce the inversion fail- 
ure. For a summary of our latest techniques, see (5lj . 



E. Diagnostics 

1. Constraints and rest-mass conservation 

During the evolution, we monitor the Hamiltonian 
and momentum constraints, calculated by Eqs. (40)-(43) 
of [H. 

When hydrodynamic matter is evolved on a fixed uni- 
form grid, our hydrodynamic scheme guarantees that the 
rest mass Mq is conserved to machine roundoff error. 
This strict conservation is no longer maintained in an 
AMR grid, where spatial and temporal prolongation is 
performed at the refinement boundaries. Hence, we also 
monitor the total rest mass, 



Mn = 



(26) 
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FIG. 3. Case A orbital-plane K contours at selected times. Contours are plotted according to K = -Km ax 10 028j ) (j—0, 1, ... 
8). Here K maK = 1.6. The color coding is the same as used in Figs. [TJand[2] A density cutoff of 10 _1 po,max has been imposed, 
where ttpo.ma* = 0.0917, or po.max = 4.58 x 10 14 g cm~ 3 (1.45M Q /M ) 2 . The dual cold core nature of the HMNS is visible and it 
becomes clear that between the two cores a hot area has formed, where 40-50% of the total pressure is due to thermal pressure. 
In the outer parts of the HMNS the contribution of the thermal component is greater than 50% of the total pressure. Here 
M = 1.32 x lO" 5 (M o /1.45M )s= 3.98(M o /1.45M )km is the ADM mass, and Mo denotes the rest mass of each star. 




1000 2000 3000 4000 1000 2000 3000 4000 O.Sfi 0.9 0.92 

t/M t/M J.nt/M 2 



FIG. 4. Left panel: Maximum rest-mass density /9o,max(t) normalized by po,max(t = 0). Middle panel: Minimum value of the 
lapse function vs time. Right panel: Maximum rest-mass density vs total angular momentum for the different cooling time 
scales considered. Here M = 1.32 x 10~ 5 (Af /1.45M Q )s= 3.98(M /1.45M Q )km is the ADM mass, and ftp ,max(0) = 0.0917, or 
Po,max(0) = 4.58 x 10 14 g cm~ 3 (1.45Mo/Mo) 2 . Mo here denotes the rest mass of each star. 



during the evolution. Rest-mass conservation is also vi- 
olated whenever po spuriously drops below and is then 
reset to the atmosphere value. This usually happens only 
in the very low-density atmosphere or deep inside the BH 
horizon where high accuracy is difficult to maintain. In 
the simulations presented in this paper the violation of 
rest-mass conservation is less than 1%. 



2. Temperature 

We measure the thermal energy generated by shocks 
via the entropy parameter K = P/P co id, where P co id = 
kpq is the pressure associated with the cold EOS. The 



specific internal energy can be decomposed into a cold 
part and a thermal part: e = e CO M + £th with 

fcoid = - J Pcoidrf(l/po) = f^Po" 1 ■ (27) 

Hence the relationship between K and e t h is 

IP k r _! 



= (K - l)e cold . (28) 

For shock-heated gas, we always have K > 1 (see Ap- 
pendix B of (33j . 

To estimate the temperature of the remnant we model 
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TABLE I. Summary of cases. The second column in- 
dicates whether cooling is applied. Here M = 1.32 x 
KT 5 (Af o /1.45M )s= 3.98(M o /1.45M )km is the ADM 
mass. 



Case Name (a) 


Cooling On Cooling 


; time scale, t/M 


A 


No 


oo 


Bl 


Yes 


150.82 


B2 


Yes 


301.64 



^ The inspiral and merger calculation is part of case A. In 
cases Bl and B2 cooling is turned on at t ~ 1600M, at 
which point the HMNS remnant of case A has relaxed to a 
quasiequilibrium state. 



the specific thermal energy as 
3/cT 

eth 



2m n 



f- 



,aT 4 



(29) 



where m n is the mass of a nucleon, k is Boltzmann's con- 
stant and a is the radiation constant. The first term rep- 
resents the approximate thermal energy of the nucleons, 
and the second term accounts for the thermal energy due 
to relativistic particles. The factor / reflects the num- 
ber of species of relativistic particles that contribute to 
the thermal energy and is determined self-consistently as 
outlined in [l4| . 

Note that the value of T depends on the adopted mass 
of the initial configuration and breaks the scale invariance 
with respect to k. For this purpose we set k = 269.6km 2 , 
for which M = 2.69M Q . 



GW extraction, energy and angular momentum 
conservation 



where e ? ;,fc is the flat-space Levi-Civita tensor. As pointed 
out in [33[ , the integrals can be evaluated more accurately 
by alternative expressions that use Gauss's law [§]: 



M, 



\ 16-7T 

16tt 24tt 



167T 



(32) 



7 _ S 1 

1 



drx^{A\ + -x 3 d n K 



-x'A km d n ^ km + 8irx j S n ) 
~e zj n j> iP 6 x j A m n dX m , 



(33) 



where S is a surface surrounding the BH horizon (when 
a BH is present), V is the volume between S and the 
outer boundary, p — n il n v T lLV ', and R is the Ricci scalar 
associated with the conformal 3-metric 7^ . 

To check the violation of energy and angular momen- 
tum conservation, we monitor the quantities 

5E = I M - M lnl! (t) - AE GW (t)\/M , (34) 
5J =\J- J int (t) - AJ GW (t)\/J , (35) 

where J and M are the ADM angular momentum and 
mass of the binary, respectively, and M- m t(t) and Jj n t(i) 
are the total mass and angular momentum of the system 
at time t as calculated by Eqs. (|32l) and (|3"3"|) . Note that 
J int (0) = J, M int (0) = M, and 5^(0) = SJ(0) = at 
t = 0. The maximum violation of energy conservation in 
the simulations we present in this paper is SE = 2% and 
the maximum angular momentum violation is S J — 3.4%. 



Gravitational waves are extracted using the Newman- 
Penrose Weyl scalar ^4 at various extraction radii be- 
tween 40M and 150M. We decompose yb^ into s = —2 
spin-weighted spherical harmonics up to and including 
I = 4 modes. At each extraction radius, the retarded time 
is computed using the technique described in Sec. IIB 
of [13 to reduce the near-field effect. 

We compute the radiated energy AEqw and z- 
component of angular momentum AJqw using expres- 
sions equivalent to Eqs. (33), (39), (40) and (49) of H. 

We also monitor the mass -Mj nt and (z-component of) 
the total angular momentum Ji nt interior to the simu- 
lation domain. These quantities are defined as integrals 
over the surface of the outer boundary dV) of the com- 
putational domain: 

Mint = ^£ v (l^ 1 ~ 7^V) (30) 
Jint = ^e k zj j x* {K% ~ 8£K)dL m , (31) 



V. RESULTS 

This section presents the results from our fully rela- 
tivistic binary NS simulations with cooling. Following 
the merger and formation of a HMNS remnant, which 
were carried out without cooling, we perform three sub- 
sequent calculations. In one calculation cooling is never 
turned on (case A) . In the other two calculations cooling 
is triggered at t w 1600M, at which point the HMNS 
remnant has relaxed to a quasiequilibrium state - the 
rest-mass density has settled and is changing on a secu- 
lar (GW) time scale. We continue the simulations with 
cooling, choosing either a short or a long cooling time 
scale, corresponding to case Bl and B2, respectively. Ta- 
ble U summarizes these different cases. 

In addition, we performed all simulations at both mod- 
erate and high resolutions. The grid configurations are 
outlined in Table |TTJ The results we obtained are insensi- 
tive to resolution, and for this reason all plots that follow 
correspond to data from our high-resolution runs. 
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TABLE II. Grid configurations. Here, JVns denotes the number of grid points covering the smallest diameter of the neutron 
star initially (13.54km = SAM). "Moderate" indicates the moderate resolution runs, and "High" the high resolution ones. 

Grid Hierarchy (in units of M)' a ' Max. resolution A?ns 

Moderate (181.98, 90.488, 45.244, 22.622, 15.081, 11.311, 7.541, 5.302) M/16.97 116 
High (181.98, 90.488, 45.244, 22.622, 15.081, 11.311, 7.541, 5.043) M/21.22 144 

There are two sets of nested refinement boxes centered on each of the NSs. This column specifies the half side length of 
the refinement boxes centered on each star. 



A. Inspiral and merger 



During inspiral and merger in case A no cooling takes 
place. The evolution of the rest-mass density contours 
in the orbital plane are shown in Fig. [TJ As gravita- 
tional waves carry angular momentum away from the 
system, the orbits become tighter, and the stars be- 
come strongly tidally distorted (top row, middle panel 
of Fig. [JJ). Shortly after the second orbit the stars col- 
lide (top row, right panel of Fig. [TJ), marking the onset of 
the merger phase. Half an orbit later, the shock-heated 
stars become strongly sheared (bottom row, left panel of 
Fig.[T} and eventually merge and settle in a quasiequilib- 
rium configuration that consists of two cold cores, sep- 
arated by hot, dense material and surrounded by a hot, 
dense mantle (bottom row, middle and right panels of 
Fig. [TJ). The upper panel of Fig. [5] shows the meridional 
rest-mass density contours at the same time as the final 
panel of Fig. [TJ No mass outflows from the system are 
observed, so the HMNS has a rest mass approximately 
equal to the initial rest mass of the system. This mass 
now exists within an equatorial radius of about 20km and 
a polar radius of about 12km. 

Fig. |3] shows the evolution of the orbital-plane K = 
P/Pcoid contours for case A at merger and following 
HMNS formation. The left panel in Fig. [3] shows how 
the collision of the two stars begins to shock heat the 
matter. In the middle and right panels in Fig.[3]the total 
pressure in the HMNS is clearly greater than the cold 
pressure everywhere except inside the two cores. Notice 
the existence of a hot area between the two cores of the 
remnant. In this area K sa 1.5, indicating that the ther- 
mal pressure adds a total of 50% additional support to 
the cold pressure. In the outer layers of the remnant 
the thermal pressure provides up to ~ 80% of the total 
pressure. The lower panel of Fig. [5] shows the meridional 
K contours. Notice that K approaches 1.8 in both the 
outer HMNS layers and the hot region between the dou- 
ble core. These results demonstrate that shock heating 
has enhanced the total pressure, which, along with cen- 
trifugal forces, contributes to the support of the remnant 
against gravitational collapse. 

Following HMNS formation at about t « 500M = 
6.6(M/2.69JVf©)ms, the remnant survives for a long qua- 
sistationary epoch, during which the maximum density 
increases almost linearly with time (see left panel of 
Fig]?]). Similar behavior is reported in [2~4j when using 



t = 2509.1M 




-40 -20 20 
X (km) 



40 



FIG. 5. Case A orbital-plane temperature contours. Contours 
are plotted according to T = T max 10~ 0136j , (j=0, 1, ... 8), 
where T max = 2.57 x lO^K = 22.18MeV. The color coding 
here is the same as in Fig. [TJ A density cutoff of 10 _2 po, max 
has been imposed, where po.max is the maximum density on 
the grid. The maximum temperature is at the center of the 
HMNS remnant. The rms temperature in the remnant is T = 



6.35 x 10 10 K » 5.5MeV. Here M 
the ADM mass. 



1.32 x 10~ b s= 3.98km is 



the same T-law EOS adopted here and no cooling. 

In addition, [241 ] performed a simulation of the same 
system, but with a strict polytropic EOS (P = k/9q), in 
which shocks are artificially suppressed. In this case, 
it is found that the resulting HMNS collapsed when 
Po,max ~ 2po,max,initiai at which point t « 21ms. Ap- 
plying this same density criterion to the T-law EOS, 
they extrapolate that the HMNS would collapse at t « 
110(M/2.69M Q )ms. Our simulations show that po,max ~ 
2/Oo,max,initiai at t « 105(M/2.69M Q )ms, in good agree- 
ment with [HJ's result. In a follow-up calculation [25jj . 
the same authors demonstrate that the HMNS remnant 
collapses to a BH at t w 13O(M/2.69M )ms. 

In the polytropic EOS simulation (shocks disallowed), 
centrifugal forces provide the only source of support 
against collapse and t co \\ ~ £gw- However, in the T- 
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law EOS simulation (shocks allowed), additional pres- 
sure from thermal support is also present, so t co \\ can be 
larger. In fact, in the absence of cooling, a sufficiently 
hot remnant may never collapse. This is what [26j con- 
clude from their case M simulation, as they demonstrated 
that a hot TOV star could support in equilibrium a mass 
greater than the total mass of their merged remnant. 

Given that the T-law EOS remnant collapses about 100 
ms later than the polytropic one, it is clear that shocks 
play a key dynamical role. However, shocks, which occur 
on a hydrodynamical time scale, give rise to at least three 
effects: 

1. they heat the gas, increasing the total pressure sup- 
port; 

2. they affect the matter profile; 

3. they redistribute the angular momentum. 

A priori it is not clear which of these effects is respon- 
sible for prolonging the HMNS lifetime, and the answer 
cannot be determined by comparing simulations that do 
not allow for shocks to those that do. 

To investigate the importance of the shock-induced 
thermal pressure support against collapse, cooling the 
hot HMNS remnant is crucial. But what is the neutrino 
cooling time scale? 

Analyzing the nascent HMNS in case A, we can now es- 
timate a realistic cooling time scale directly from our sim- 
ulation data. In Fig. [5] we show the case A orbital-plane 
temperature contours of the HMNS remnant. The maxi- 
mum and rms temperatures are ~ 22MeV and ~ 5.5MeV, 
respectively. Using these values for the neutrino energy, 
and setting M = 2.7M e and R = 16km, Eq.© yields a 
neutrino diffusion time scale of ~ 165ms — 2.64s. Note 
that this range is consistent with our discussion in Sec. HTl 
where the neutrino diffusion time scale was estimated to 
be comparable to the gravitational wave time scale of 
- 120ms. 

Therefore, as the additional thermal pressure is a dy- 
namically important source of support, then cooling must 
be incorporated in numerical simulations to accurately 
determine the time interval between merger and delayed 
collapse. Low energy neutrinos (E v < lOMcV), that es- 
cape the HMNS on a time scale < tow, can remove a 
sufficiently large amount of thermal energy to accelerate 
the collapse. This is the main point we want to emphasize 
in this paper. Our conclusion may even be more impor- 
tant when a stiff EOS is employed, as the HMNS will be 
less compact resulting in less effective GW emission. 



B. Cooling study 

To assess the impact of the shock-induced thermal 
pressure in the HMNS remnant, we slowly remove the 
thermal pressure using our covariant cooling technique. 
We chose two cooling time scales for our study. The 



short one is r c i = 6.5id for case Bl and the long one is 
t c ,2 = 2r c i = 13t<i for case B2, where 

i d = 4= ~ 23M - ( 36 ) 



is the dynamical time scale of the HMNS, where p is the 
mean HMNS density. This choice is arbitrary, but it is set 
so that the cooling time scale is significantly longer than 
the dynamical time scale of the remnant, as in physically 
realistic stars, but short enough for our simulations to be 
completed within a reasonable time. 

Our runs with cooling lead to BH formation within a 
few cooling time scales. This can be seen in the left two 
panels of Fig. U where we plot the maximum rest-mass 
density po.max and the minimum of the lapse function 
own versus time, respectively. From these plots it be- 
comes clear that case A does not collapse within the in- 
tegration time, while both cases Bl and B2 form a BH. 
Notice that in cases Bl and B2, when po.max roughly 
equals two times its initial value, the lapse function col- 
lapses and a BH forms. Note that this is consistent with 
the polytropic runs of 24], in which shocks were sup- 
pressed. Therefore, in all cases with the adopted EOS, 
collapse takes place when p ,max ~ 2p ,max,initiai, when 
thermal energy is drained from the system. 

If the collapse is driven by cooling, then we naturally 
expect that a longer cooling time scale will increase the 
lifetime of the HMNS. This is precisely what we find: the 
collapse in case B2 occurs later than in case Bl. Further 
evidence of cooling-induced collapse is shown in the right 
panel of Fig. 0] There the tracks of po.max against the 
total angular momentum of the remnant J- m t are plot- 
ted using Eq. (|33p . This plot demonstrates that during 
the post-merger evolution for the same Ji n t, po,max is al- 
ways larger with cooling present (e.g. po,max 

in cases Bl, 

B2 is ~ 30% larger than in case A for the smallest J; n t 
reached in case A). Therefore, it is the reduction of ther- 
mal pressure that leads to a more compact remnant and 
not angular momentum loss driven by GWs. 

In Fig. [HI we show the evolution of the rest-mass den- 
sity (upper row) and K contours (lower row) for case B2. 
The selected times correspond to 1, 8, and 11 cooling 
time scales after cooling was turned on (54[. The right 
panel corresponds to a time shortly before an apparent 
horizon forms. These plots indicate that as thermal pres- 
sure is removed, the double cores approach one another 
and the HMNS becomes more compact (notice the den- 
sity increase and the shrinking size of the remnant with 
increasing time). We find that when the hot area between 
the two cores is cooled to if « 1.05 the two cores merge 
and form a single-core HMNS. Shortly after this occurs 
(at t « 3000Af) the remnant undergoes catastrophic col- 
lapse. 

Based on these results we conclude that thermal pres- 
sure contributes significantly to support against collapse. 
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t=1802.3M t=2509.1M t=3270.3M 
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FIG. 6. Upper row: Case B2 orbital-plane rest-mass density contours at selected times. Contours are plotted according to 
po = po,max(10-°- 2j+0 ' 568 ), (j=0, 1, ... 8), where K p , m ax = 0.0917, or p , ma x = 4.58 x 10 14 g cm- 3 (1.45M Q /M ) 2 . Lower row: 
Case B2 orbital-plane K contours at selected times. Contours are plotted according to K = -ftT m axlO _0 ' 025j , (j=0, 1, ... 8). 
Here ]f max = 1.6. The color coding is the same as used in Fig. [T] with light blue indicating to K w 1, yellow K ~ 1.2 and dark 
red K sa 1.4. Here M = 1.32 x lO" 5 (M o /1.45M )s= 3.98(M /1.45M©)km is the ADM mass, and M denotes the rest mass of 
each star. 



C. Angular momentum conservation 

Figure [7] plots the evolution of total angular momen- 
tum (Ji n t) and angular momentum carried off by GWs 
(AJgw)j normalized by the ADM angular momentum 
of the binary (J) for cases A and B2. We find that 
(Jint + AJgw)/<A which should be equal to unity at all 
times, is the same and close to unity for all three cases A, 
Bl and B2. This implies that cooling carries off negligible 
amounts of angular momentum which is consistent with 
earlier estimates |55| . The maximum violation of angular 
momentum conservation is SJ rs 3.4%. Notice that Ji„t 
is smaller in case B2 than in case A, while A Jew is larger 
in case B2 than in case A. The distinction is due to the 
fact that as thermal energy is radiated away the remnant 
becomes more compact, enabling GWs to remove angular 
momentum faster. We conclude that cooling accelerates 
the collapse of the HMNS by the combined action of two 
effects: 

1. Cooling removes thermal pressure support, yielding 



a more compact remnant. 

2. As the remnant becomes more compact, GWs are 
able to carry away angular momentum more effi- 
ciently. 

VI. SUMMARY AND FUTURE WORK 

A differentially rotating, quasiequilibrium HMNS is 
a transient configuration that can arise following the 
merger of NSNS binaries. The mass of a HMNS is larger 
than the maximum mass that can be supported by a cold 
EOS, even with maximal uniform rotation. A HMNS will 
eventually undergo "delayed collapse" on a secular (dis- 
sipative) time scale and may power a sGRB. 

When HMNSs are born in NSNS mergers, they are 
rapidly differentially rotating and hot due to shock heat- 
ing. Therefore, HMNSs will collapse to a BH either on an 
angular momentum loss/magnetic braking time scale or 
on a cooling time scale. A priori it is not clear which of 
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FIG. 7. Angular momentum vs time. Here Jmt(A) and 
AJgw(^4) are the total angular momentum and the angu- 
lar momentum carried away by GWs, respectively, for case 
A. Jint(-B2) and AJgw(S2) correspond to case B2. The sum 
in both cases is the same and given by Jmt + AJgw- All 
quantities are normalized by the ADM angular momentum of 
the binary J. 



the two above mechanisms is most important for holding 
up an HMNS against collapse: centrifugal forces or ther- 
mal pressure. The answer to this question is still open 
and may depend on the stellar model, EOS, and initial 
magnetic fields. 

Determining which mechanism drives a HMNS to col- 
lapse has observational consequences; the time scale of 
collapse will set the interval between the NSNS merger 
chirp signal and the delayed collapse burst signal, which 
may be measured by LIGO/VIRGO. Careful modeling 
of HMNS physics will thus place constraints on magnetic 
field magnitudes, the existence of bar modes, and/or the 
relevant cooling mechanisms. In addition, such observa- 
tions could place constraints on the temperature of mat- 
ter as well as the nuclear EOS. 

To disentangle the effects of thermal support from 
those of rotational support, previous studies compared 
results from NSNS simulations that suppress shocks to 
those that allow shocks. If the HMNS remnant lives 
longer in the case with shocks than without, then it is 



tempting to infer that thermal pressure due to shock 
heating is solely responsible. However it is not possi- 
ble to draw such a firm conclusion because shocks, which 
act on a hydrodynamical time scale, not only heat the 
gas, thereby increasing the total pressure support, but 
also affect the matter profile and redistribute angular 
momentum. Different matter and angular momentum 
profiles alone can increase the lifetime of a HMNS via an 
increase of both the GW time scale and the amount of 
differential rotational support. 

To address this issue, we first performed long-term, 
high-resolution GRHD NSNS simulations through inspi- 
ral, merger, and HMNS formation, allowing for shocks. 
Following HMNS formation, we continue the evolution 
both with and without cooling. When cooling is turned 
off, the remnant collapses on the GW time scale. How- 
ever, when cooling is turned on we find that the HMNS 
collapses and forms a BH within a few cooling time scales. 

Our simulations demonstrate that shock-induced ther- 
mal pressure is a significant source of support against 
gravitational collapse in the case of a stiff L-law EOS - 
a result consistent with simulations that employ a more 
realistic EOS [26| - and show explicitly that cooling can 
induce the catastrophic collapse of a HMNS. Estimating 
the temperature of the HMNS remnant, we find that a re- 
alistic neutrino cooling time scale is of order a few 100ms. 
Given that the estimated cooling and angular momen- 
tum loss/magnetic braking time scales can be compara- 
ble, cooling should be accounted for to accurately deter- 
mine the lifetime of a HMNS. Therefore simulations that 
implement cooling will lead to earlier collapse than sim- 
ulations that ignore it, otherwise the predicted GW and 
EM signatures from these delayed collapse events may be 
incorrect. 

Therefore, to accurately determine the lifetime of 
HMNS remnants, neutrino cooling physics should be in- 
corporated in NSNS simulations. In the future we plan to 
revisit the subject using a more realistic neutrino leakage 
scheme, such as that used in (26l.[56j. in conjunction with 
more realistic treatment of the microphysics involved. 
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